#include <stdio.h>
#include <math.h>
#include "eno.h"

#define data_max_x 132

#define pi 3.14159265358

double z1[data_max_x+1][2*data_max_x+1],
       z2[data_max_x+1][2*data_max_x+1];
int num_side1,num_side2,i,xstep,ystep,boundary_size;
FILE *in_file1,*in_file2;
char in_file_name1[20],in_file_name2[20];
double z_avg,l1_norm,l2_norm,temp_float,r_end,hh;

int slide_number1,slide_number2,num_side1_y,num_side2_y,use_distance_int,
    slide_count;
double file_var,density_above,density_below,current_delta_x,all_points;
logical use_distance;


double smooth_front(in_distance,value_above,value_below,
                 value_spread,smooth_outside)
double in_distance,value_above,value_below,value_spread;
logical smooth_outside;
{
double cutoff,slope,intercept,range;

  cutoff = value_spread*current_delta_x;
  if (smooth_outside)
    if (value_below > value_above)
      in_distance -= cutoff;
    else
      in_distance += cutoff;

  intercept = (value_above+value_below)/2.0;
  range = (value_above-value_below)/2.0;

  if (cutoff > 0.0)
    slope = range/cutoff;

  if (in_distance == 0.0)
    return intercept;
  else if (in_distance < -cutoff)
    return value_below;
  else if (in_distance < cutoff)
/*  return (slope*in_distance + intercept); linear smoothing */
    return (intercept + range*sin(pi*in_distance/(2.0*cutoff)) );
  else if (in_distance > cutoff)
    return value_above;
}

double get_r(in_xstep,in_num_side)
int in_xstep,in_num_side;
{
  if (use_distance_int==2)
   return r_end*(in_xstep-0.5)/in_num_side;
  else
   return 1.0;
}

double cvt_density(in_distance)
double in_distance;
{
double density_spread;

  density_spread=1.5;

  if (use_distance)
    return smooth_front(in_distance,density_above,density_below,
            density_spread,FALSE); 
  else
    return in_distance;
}

void main()
{
int scan_status;

  printf("This program gives l2 norm of the difference between\n");
  printf("two files where the finer grid is staggered.\n");
  printf("contents of file : u at x1,y1 ... xn,y1\n");
  printf("                        x1,y2 ... xn,y2 ...\n");
  printf("First entry in file is num_side_x.\n");

  printf("Enter boundary size (e.g. 0,1,2) =>\n");
  scan_status = scanf("%d",&boundary_size);
 
  printf("Enter file 1 (coarse) => \n");
  scan_status = scanf("%19s",in_file_name1);
  printf("Enter slide number1 (0..n) =>\n");
  scan_status = scanf("%d",&slide_number1);

  printf("Enter file 2 (fine) => \n");
  scan_status = scanf("%19s",in_file_name2);
  printf("Enter slide number2 (0..n) =>\n");
  scan_status = scanf("%d",&slide_number2);

  printf("Enter:  1-distance function   2-axisymmetric=>\n");
  scan_status = scanf("%d",&use_distance_int);
  use_distance = (use_distance_int==1);
  if (use_distance_int==2) 
   {
    printf("Enter r_end=>\n");
    scan_status=scanf("%lf",&r_end);
    printf("r_end is %10.5f\n",r_end);
   }

  if (use_distance)
    {
     printf("enter density_above (outside bubble)=>\n");
     scan_status=scanf("%lf",&density_above);
     if (density_above<500.0) 
       density_below=1000.0;
     else 
       density_below=1.0;
    }
  else
    {
     density_above=density_below=1.0;
    }
  printf("density_above,density_below%10.4f%10.4f\n",
    density_above,density_below);

  in_file1 = fopen(in_file_name1,"r");
  scan_status = fscanf(in_file1,"%d\n",&num_side1);
  printf("coarse grid : %5d\n",num_side1);
  printf("enter number of points in y direction for coarse grid=>\n");
  scan_status = scanf("%d",&num_side1_y);

  for (slide_count=0;slide_count<=slide_number1;slide_count++)
   {
    current_delta_x=7.0/num_side1;
    for (ystep = 1;ystep <= num_side1_y;ystep++)
     for (xstep = 1;xstep <= num_side1;xstep++)
      {
       scan_status = fscanf(in_file1,"%lf\n",&file_var);
       z1[xstep][ystep]=cvt_density(file_var);
      }

    if (slide_count<slide_number1)
      scan_status = fscanf(in_file1,"%d\n",&num_side1);
    printf("it side %5d\n",num_side1);
   }
  scan_status = fclose(in_file1);
  
  in_file2 = fopen(in_file_name2,"r");
  scan_status = fscanf(in_file2,"%d\n",&num_side2);
  printf("fine grid : %5d\n",num_side2);
  num_side2_y=num_side1_y*(num_side2/num_side1);

  for (slide_count=0;slide_count<=slide_number2;slide_count++)
   {
    current_delta_x=7.0/num_side2;
    for (ystep = 1;ystep <= num_side2_y;ystep++)
     for (xstep = 1;xstep <= num_side2;xstep++)
      {
       scan_status = fscanf(in_file2,"%lf\n",&file_var);
       z2[xstep][ystep]=cvt_density(file_var);
      }
    if (slide_count<slide_number2)
      scan_status = fscanf(in_file2,"%d\n",&num_side2);
    printf("it side %5d\n",num_side2);
   }
  scan_status = fclose(in_file2);

  l2_norm = 0.0;
  l1_norm = 0.0;
  for (ystep = 1+boundary_size;ystep <= num_side1_y-boundary_size;ystep++)
    for (xstep = 1+boundary_size;xstep <= num_side1-boundary_size;xstep++)
      {
       if (num_side2 > num_side1)
         {
          z_avg = z2[2*xstep-1][2*ystep-1]+z2[2*xstep][2*ystep-1]+
                  z2[2*xstep][2*ystep]+z2[2*xstep-1][2*ystep];
          z_avg /= 4.0;
         }
       else
         z_avg = z2[xstep][ystep];
       z_avg = fabs(z_avg-z1[xstep][ystep]);
       if (use_distance_int==2)
        z_avg *= 2.0*pi*get_r(xstep,num_side1)*r_end/num_side1;
       l2_norm += z_avg*z_avg;
       l1_norm += fabs(z_avg);
      }  

  l2_norm = sqrt(l2_norm);

  if (use_distance_int==2)
   {
    l1_norm/=num_side1_y*pi*sqr(1.5);
    l2_norm/=sqrt(num_side1_y*pi*sqr(1.5));
   }
  else
   {
    all_points=(num_side1-2*boundary_size)*(num_side1_y-2*boundary_size);
    l1_norm/=all_points;
    l2_norm/=sqrt(all_points);
   }
  printf("l1_norm,l2_norm %12.8f%12.8f\n",l1_norm,l2_norm);
  return;
}
